library(mi)
library(nnet)
library(foreign)
library(Amelia)
library(norm)
library(mi)
library(nnet)

source("realhelper.R")

part <- 1
resfile <- paste(c("results",part,".csv"), collapse="")

results <- as.numeric()
resnames <- c("iteration", "method", "imputation_method", "time", 
	"mvmatch.co", "mvmatch.bias.co", "coef.co", "coef.bias.co", "coef.mns.co", "fit.all.co", "fit.all.bias.co", "fit.av.co", "fit.av.bias.co", 
	"mvmatch.bi", "pr.rmse.bi", "pr.bias.bi", "coef.bi", "coef.bias.bi", "coef.mns.bi", "fit.all.bi", "fit.all.bias.bi", "fit.av.bi", "fit.av.bias.bi", 
	"mvmatch.or", "pr.rmse.or","pr.bias.or", "coef.or", "coef.bias.or", "coef.mns.or", "fit.all.or", "fit.all.bias.or", "fit.av.or", "fit.av.bias.or", 
	"mvmatch.un", "pr.rmse.un","pr.bias.un", "coef.un", "coef.bias.un", "coef.mns.un", "fit.all.un", "fit.all.bias.un", "fit.av.un", "fit.av.bias.un")

for(it in 1:20){
	iter <- 20*(part-1) + it
	for(im in c("ppd", "pmm")){
		dgp <- .dgp(miss.prop=0.1, imp_meth=im, use.miss=FALSE)
		
		res <- .amelia(dgp, m=5, command="amelia")
		results <- rbind(results, c(iter, "Amelia", im, res$time,
			mean(res$co$mvmatch), mean(res$co$mvmatch.bias), res$co$coef, res$co$coef.bias, res$co$coef.mns, 
			res$co$fit.all, res$co$fit.all.bias, res$co$fit.av, res$co$fit.av.bias, 
			mean(res$bi$mvmatch), mean(res$bi$pr.rmse), mean(res$bi$pr.bias), res$bi$coef, res$bi$coef.bias, res$bi$coef.mns, 
			res$bi$fit.all, res$bi$fit.all.bias, res$bi$fit.av, res$bi$fit.av.bias, 
			mean(res$or$mvmatch), mean(res$or$pr.rmse), mean(res$or$pr.bias), res$or$coef,res$or$coef.bias, res$or$coef.mns, 
			res$or$fit.all, res$or$fit.all.bias, res$or$fit.av, res$or$fit.av.bias, 
			mean(res$un$mvmatch), mean(res$un$pr.rmse), mean(res$un$pr.bias),res$un$coef, res$un$coef.bias, res$un$coef.mns, 
			res$un$fit.all, res$un$fit.all.bias, res$un$fit.av, res$un$fit.av.bias))
		
		res <- .amelia(dgp, m=5, command="norm")
		results <- rbind(results, c(iter, "Norm", im, res$time,
			mean(res$co$mvmatch), mean(res$co$mvmatch.bias), res$co$coef, res$co$coef.bias, res$co$coef.mns, 
			res$co$fit.all, res$co$fit.all.bias, res$co$fit.av, res$co$fit.av.bias, 
			mean(res$bi$mvmatch), mean(res$bi$pr.rmse), mean(res$bi$pr.bias), res$bi$coef, res$bi$coef.bias, res$bi$coef.mns, 
			res$bi$fit.all, res$bi$fit.all.bias, res$bi$fit.av, res$bi$fit.av.bias, 
			mean(res$or$mvmatch), mean(res$or$pr.rmse), mean(res$or$pr.bias), res$or$coef,res$or$coef.bias, res$or$coef.mns, 
			res$or$fit.all, res$or$fit.all.bias, res$or$fit.av, res$or$fit.av.bias, 
			mean(res$un$mvmatch), mean(res$un$pr.rmse), mean(res$un$pr.bias),res$un$coef, res$un$coef.bias, res$un$coef.mns, 
			res$un$fit.all, res$un$fit.all.bias, res$un$fit.av, res$un$fit.av.bias))
		
		res <- .completecase(dgp)
		results <- rbind(results, c(iter, "Complete Cases", im, NA,
			NA, NA, res$co$coef, res$co$coef.bias, res$co$coef.mns, 
			NA, NA, res$co$fit.av, res$co$fit.av.bias, 
			NA, NA, NA, res$bi$coef, res$bi$coef.bias, res$bi$coef.mns, 
			NA, NA, res$bi$fit.av, res$bi$fit.av.bias, 
			NA, NA, NA, res$or$coef,res$or$coef.bias, res$or$coef.mns, 
			NA, NA, res$or$fit.av, res$or$fit.av.bias, 
			NA, NA, NA,res$un$coef, res$un$coef.bias, res$un$coef.mns, 
			NA, NA, res$un$fit.av, res$un$fit.av.bias))
		
		res <- .mi(dgp, m=5, n.iter=0, n.chains=5, est="MNL")
		results <- rbind(results, c(iter, "MCAR", im, res$time,
			res$co$mvmatch, res$co$mvmatch.bias, res$co$coef, res$co$coef.bias, res$co$coef.mns, 
			res$co$fit.all, res$co$fit.all.bias, res$co$fit.av, res$co$fit.av.bias, 
			res$bi$mvmatch, res$bi$pr.rmse, res$bi$pr.bias, res$bi$coef, res$bi$coef.bias, res$bi$coef.mns, 
			res$bi$fit.all, res$bi$fit.all.bias, res$bi$fit.av, res$bi$fit.av.bias, 
			res$or$mvmatch, res$or$pr.rmse, res$or$pr.bias, res$or$coef,res$or$coef.bias, res$or$coef.mns, 
			res$or$fit.all, res$or$fit.all.bias, res$or$fit.av, res$or$fit.av.bias, 
			res$un$mvmatch, res$un$pr.rmse, res$un$pr.bias,res$un$coef, res$un$coef.bias, res$un$coef.mns, 
			res$un$fit.all, res$un$fit.all.bias, res$un$fit.av, res$un$fit.av.bias))
			
		res <- .mi(dgp, m=5, n.iter=30, n.chains=5, est="MNL")
		results <- rbind(results, c(iter, "MI: MNL", im, res$time,
			res$co$mvmatch, res$co$mvmatch.bias, res$co$coef, res$co$coef.bias, res$co$coef.mns, 
			res$co$fit.all, res$co$fit.all.bias, res$co$fit.av, res$co$fit.av.bias, 
			res$bi$mvmatch, res$bi$pr.rmse, res$bi$pr.bias, res$bi$coef, res$bi$coef.bias, res$bi$coef.mns, 
			res$bi$fit.all, res$bi$fit.all.bias, res$bi$fit.av, res$bi$fit.av.bias, 
			res$or$mvmatch, res$or$pr.rmse, res$or$pr.bias, res$or$coef,res$or$coef.bias, res$or$coef.mns, 
			res$or$fit.all, res$or$fit.all.bias, res$or$fit.av, res$or$fit.av.bias, 
			res$un$mvmatch, res$un$pr.rmse, res$un$pr.bias,res$un$coef, res$un$coef.bias, res$un$coef.mns, 
			res$un$fit.all, res$un$fit.all.bias, res$un$fit.av, res$un$fit.av.bias))
			
		res <- .mi(dgp, m=5, n.iter=30, n.chains=5, est="RNL")
		results <- rbind(results, c(iter, "MI: RNL", im, res$time,
			res$co$mvmatch, res$co$mvmatch.bias, res$co$coef, res$co$coef.bias, res$co$coef.mns, 
			res$co$fit.all, res$co$fit.all.bias, res$co$fit.av, res$co$fit.av.bias, 
			res$bi$mvmatch, res$bi$pr.rmse, res$bi$pr.bias, res$bi$coef, res$bi$coef.bias, res$bi$coef.mns, 
			res$bi$fit.all, res$bi$fit.all.bias, res$bi$fit.av, res$bi$fit.av.bias, 
			res$or$mvmatch, res$or$pr.rmse, res$or$pr.bias, res$or$coef,res$or$coef.bias, res$or$coef.mns, 
			res$or$fit.all, res$or$fit.all.bias, res$or$fit.av, res$or$fit.av.bias, 
			res$un$mvmatch, res$un$pr.rmse, res$un$pr.bias,res$un$coef, res$un$coef.bias, res$un$coef.mns, 
			res$un$fit.all, res$un$fit.all.bias, res$un$fit.av, res$un$fit.av.bias))	
	}
}

colnames(results) <- resnames
write.table(results, file=resfile, sep=",", col.names=TRUE, row.names=FALSE)